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ABSTRACT 



The equations of motion are developed and solved numerically 
for the trajectory of a spherical particle passing tlirough a Prandtl- 
Meyer ejq^ansion fan. The effect of a change in y and n is shown, 
where y is the ratio of specific heats of the gas, and n/2 is the 
e>q3onent in 

, n/2 

= ( T^- ) 

o 

the assumed viscosity temperature relationship. 

It is demonstrated that for particles the order of a micron in 
diameter, slip flow will exist, and a proposed correction to the drag 
to account for this discontinuous nature of the flow is investigated. 

The results are plotted showing particle trajectory profiles 
and the components of the relative velocity or slippage velocity of 
the particle for flow deflection angles up to 70 degrees. 
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I. INTRODUCTION 

.Heterogeneous flow of a solid particle - gas mixture has been 
a topic of increased interest in recent years. The addition of light 
metals to solid roclcet fuels has focused attention on the non- equilibri- 
um nature of two-phase flovt'S in expansion processes, and numerous 
papers have been written regarding this phenomencni. 

Sarly studies of t%vo-phase flows were involved with investi- 
gations of v/ater droplet trajectories around a variety of body shapes 
in connection v/ith icing of aircraft surfaces. More recently, hetero- 
geneous nozzle flow has been examined. Gilbert, Davis, and Altman 
(reference i) related roclcet thrust losses to solid particle sizes and 
solved linearized nozzle equations assuming Stokes drag law. Kliegel 
(reference Z) treated simultaneous velocity and thermal lag along witli 
the effect of particle lag on gas properties in a one-dimensional noz- 
zle. Rannie (reference 3) showed that for moot cases of interest, 
results of equal accuracy to niimerical integration can be obtained 
through a perturbation analysis of temperature and velocity lags in a 
one -dimensional nozzle. Rannie also introduced slip flow corrections 
for drag and heat transfer. An excellent stimmary of recent advances 
in heterogeneous nozzle flow was presented at the American Rocket 
Society, Solid S-ropcUant Rocket Conference in January 1962, by S.. 

F. Hoglund (reference 4). 

Several problems relating to supersonic flow of tv^o-phase 
systems have also been examined. Carrier (reference 5) investigated 
a solution to particle flov/ tlirough normal shocks, and oblique shocks 
were treated by Morgenthaler (reference 6). 
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The behavior of a single particle passing through a Prandtl- 
Meyer ejqjansion fvirnishes an e?:a.mple of heterogeneous flow that can 
be analysed in a straightforward manner, and serves to illustrate 
some of the problems and important parameters of two-phase flow in 
general. This paper presents the development of the equations of mo- 
tion of small spherical particles flowing with a gas that ejq>cricnces a 
Prandtl-Meyer ejq>ansion, and tlie numerical solution of these equa- 
tions. The effects of the temperature dependence of viscosity, cloanges 
in the ratio of gas specific heats, y , and slip flow correction to the 
drag coefficient v/cre investigated by varying these factors in succes- 
sive solutions. Tlie results of this investigation have direct bearing 
on the exiiaust of aii undcrexpanded rocket iiozzle as the flow expands 
around the lip of the nozzle. 
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II. DEVELOPMENT OF BASIC ECUATIONS 

An ideal two-dimensional prandtl- Meyer expansion about a 
sharp corner in a semi-infinite channel is the assumed model for 
this investigation. A cylindrical coordinate system with origin at 
the vertex of the expansion fan is used to describe the position of a 
spherical particle passing through the fan. Denoting the radial dis- 
tance from the origin to the center of the particle as , and the 
angxilar displacement of the particle from the initial characteristic 
as (Fig- 1)» the conservation of momentum equations can be 
written 

r. QllrfO!' h -f 6 7rcr^ (U 3 - r) (i) 



[m r- O) 

dt ^ 



■f (oVr (T - re)r 



( 2 ) 



It is assumed that the only force acting on the particle is a drag force 
which can be represented by Stokes drag law modified by the factor 
^ . The radial and tangential components of the gas velocity are 
and 1/^ respectively, and is the local value of the gas viscosity. 
The particle has mass m and radius cr . 

It is convenient to rewrite the equations with ©■ the independ- 
ent variable. Choosing L - m j the angular momentum, and 
X- as the desired dependent variables, the time derivative can be 
written 




) 



) . 



Equation 1 then becomes 



f GTT u ^ 



\ — 


m 1 




d r ■ 




'de , 


Imr^J 




m r ^ 






Putting 



and mxiltiplying by 



dr ^ d('/r) 

de d e 

d'C/r) , , d('/r) ^dL 



m r' 



1 ^ d 'I /r; . I d r/r; d L \ 

^ 6 9^ 6^ ^6b' 



+ 



r 






p f 6 vra-ju mr^Uj f ^ Tr<r/J. L _ q 

cf O 



Equation 2 is simply 

dL z' m 'xT^ 



= f 6TTcr/^ ("- 



- r‘ 



J . 



(3) 



(4) 



d & ' ' ' ^ L 

Replacing d L in. the second term of equation 3 by equation 4 we 



have 



d O 

2 d^(/r) 






-f t>- r 6Tr(T mr"’A/'o 

r c/e 



c/ ("/^) 



de 



b f GTTcr /?7 ^ 



(5) 



Non-dimensional parameters are defined as 

r L 



R = 



5 



and 



A = 



( 6 ) 



m c A 

where X ~ ~ , y is the ratio of gas specific heats, and C 

is the theoreticeil maximum expeinsion velocity determined by the 
equation 



- 2CpT, ( I p 



x'-i 



m/) 



f> -O I ■ ■ 2 

The characteristic length, ^ , has the value 



(7) 






m c 

(dTr<r/j, 



(8) 




i 
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The subscript zero indicates initial conditions of the gas or particle 
as they enter the expansion fan. 

Introducing equations 6 into equation 5, we obtain 






^Ku 

d& 3 



= O 



and dividing by y ^ iS TT <J~ rn yU c 



C , 



A" 



Rj 



i- 



(_ C of ^ c J 



Equation 4 can be written in terms of the dimensionless parameters 
as 

=f f: 



^ ^ _ ) ) 
AA c / 



^ = -T/?^ ^ r A? \ ] 



( 10 ) 



The conservation equations of frictionless flow yield , 

c 

.y^ , and — as functions of S- . The equations may be writ- 

C 77 

ten in terms of our assumed coordinate system as 



Continuity: ^ ^ -t J t* 



da 



0 



Momentum: 



15. _ 35. = o 

r d& r 



'Vq 



r 



LL, + 






(H) 

( 12 ) 

(13) 



del d e 



(14) 



Energy: ^ q 2 . _ ^ ^ 

Equations 11 and 13 can be recombined to eliminate the term in the 
parentheses and show that 

Vg — a . (15) 

Using this resxilt in equation 14, Vj can be found as a fxinction of 
^3 which can, in turn, be used in equation 12 to give 

^^-=5in(A0+-o() ^ 

( 16 ) 



c 



— =: XCOs(X6+-'>^) . 

C ^ ^ 

The value of the constant oi can be found by evaluating equations 16 
at the initial conditions O- — O ^ O.^ ^ and using equation 

14 to give 

1 

2 






(17) 



The quantity 



can be found by applying equation 14 to show that 

• o 



To 



'xn- (^): 



or by equations 16 






t; 



Cos C< 



(18) 



It is assumed that viscosity varies with temperature raised to a frac- 
tional power, hence 
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or from equation 18 



it. — CO£) ^ ( X <^ 



Ut 



cos’'" o< 



(20) 



where n will normally be between 1. 0 and 2. 0. 

The equatioiEof particle motion (9 and 10) can nov/ be written 
entirely in terms of the dimensionless variables R , , the 

drag correction factor , and the constants <^< and A . Letting 



( )^ denote £1 



j ( ) 



d t> 



they become 



A' 



f' I \ " 



Al-fe) 



+- 



, ^2 ■:A.."iA0ro()r^ ^ ^ 

-f 1 - [ r-. ; - -■ ■ , /• ) 



•f 



f.\C' T-U) 



(21) 



A ' P ' 



c - s^(.\ 0 +-<x) 






COo^'cX 



R 1 

A COS(A9- - Aj =. O , (22) 



la 



A value of jl for spheres in the subsonic slip flow regime has 
been proposed by V. . D. Rannie (ref. 3 ) which, for relative Mach num- 
bers less thaii unity and for Reynolds numbers less than about 0. 5, 
taltes the form 



-;r - 1 . 'T) -H 0 £ ./ ( 5 ' ’ 

}- ^ A'i* 



A / 



(23) 



where A//^ and Rq are the relative >'^ach number and Reynolds 



number respectively. 

■ 

Re 

path to the sphere diameter, can be evaluated as follows; 



, which is related to the ratio of moleciilar mean free 



8' 



(24) 



S± ^ A- fJL](^](A] 



FroiYi equations 15 and 16 

Qo _ Coso<- 
CL ■“ cos(xe + c<) 



(25) 



Isentropic relationship betv/een density and temperature along with 



equation 19 gives 



P 



COS 

COS(X© + °<} 



_ 2 _ 

Y-1 



(26) 



Using equations 20, 25, and 26 in equation 24, we get 



Re 



a»cr 



COScK 

cos (X ^ + °< ) 



x^ 



n 

1 



(27) 



and equation 23 becomes 



■f 



I + 



z.ozfT A 

Po 



C05°< ■' 

COS ( X © + 0() 



1 



- n 



(28) 



The drag correction factor cannot be written in terms of the dimen- 
sionless similarity parameters R. and A , but is dependent on the 
particle size and initial conditions. 

It is convenient to relate the trajectory of the particle to that 
of a particle having zero relative velocity, i. e. , moving along a gas 
streamline. Designating this path as ^ ^ ^ write 



- U - Al 

di ^ d O d t ’ 
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and since 



de _ 




dt 




d r, _ 


r, u 


de 


^3 



v/e have 



Integrating this equation we obtain 

R< 






^9 
Ro 



cos oc 



_l_ 

A’- 



R. \cos(Ae 4 -o<)y 
From the definitions of and A it follows that 



Vo R’o 



A, 



which, when combined with equations 16 and 30, gives 



A 5 _ / COS 

Ao ~ lcos(xe+ o<) 



( 29 ) 



(30) 



(31) 



Two new dimensionless parameters are defined 




from which it follows that 



R 



2^ 

Ro 



cos ('A © + 
c os cx 






% 



(32) 



(33) 






/ 

(-b) ^ - 


' COSiXB -t-oOlA* 


x'- ^ tan(x & -hx) 


VR/ R, 


COSX J 





( 34 ) 
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I 



/_Lf _ _L ' 


' cos (X® + «x)' 


x^’ 


\R/ Re 


Cos oi. 


L- 



2r 




tan (\0 + o( ) 



-¥ 




^(X9^o() - ij 



( 35 ) 



A 



YAo 



COS o( 

c 05 ( xe + o(). 




(36) 



A = A, 



CC>5 o< 



cos(x e + oc)_ 






y + y(-^ - X)tan(xe-t-o() 



(37) 



AJiso from equations 6 and 16 
Ao = Rc COS - 



(38) 



Writing equations 21 and 22 in terms oi X , Y parameters 
using equations 33 through 38, we obtain the result 

X an^^A 0 ~ "t an + «X.) 

fx' 

^ ^ '^^“U^yv(c o 5(xe ^ ' 



( 39 ) 
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y (y - >') tan(Ae+o<) 



-f XR .(- 



-X y 

x^y 



f c O Sc^ 

rco5(Xe + o()U’- ^40) 



It is assumed that the particle is moving with no relative velocity prior 
to reaching the first characteristic line of the expansion. The initial 

conditions are then 

X (0 ) zr i . O , 

V (O) - 1 O , 

X'(0) = O • 

Equations 28, 39, 40, and 41 constitute the complete statement of 
the problem. 




ii 
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m. AUXILIARY EQUATIONS 

In order to specify the trajectory, the velocity as well as the 
position of the particle should be luiovvn. The important parameters 
in terms of X and Y are then 

^ = A . E 5 = X y , (42) 

R Ag 



^ - B. .El ~ I r COS • 
Rc ~ R’s Ro x[cosCA6+o() 






(43) 



dR 



A _ A - A ' 

^3 Lig dt V3 U5 R de 



which, by 33, 34, and 16, is 



^ = xy 



— — X c o"t ( A 0 + ) 

X 



(44) 



The relative velocity or slip velocity of the particle relative 
to the absolute gas velocity is simply 



VI - -y 

^3 

ag - u 

U g L( g 

The relative Mach number is then 



^3 



( 45 ) 



Ms = / / "t" 



CL 



Ms — 




(46) 
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IVo SCL-UTION OF EOUATIOMS 

Since equations 39 and 40 are non-linear, a numerical solu- 
tion is indicated for arbitrary values of Ro • For cither * 

or I the possibility of a small perturbation solution 

e:dsts; however, this possibility was not investigated. 

Equations 21 and 22 were first solved utilizing a desk calcu- 
lator and a simple iteration procedure where 




A "O' 






(47) 



— Ai 4- 4- 



L (a:) 



A Re 



\A-=- A, 



•t- A ©') 



A© 



For tliese calculations it was assumed that -i. o (Stokes drag lav/ 
applies), ■S'-=.l.4 , n - 20 , and Mo - • Various values 

of Ro between 0. 1 and 100 w^ere investigcited; however, it was 
found that the iteration interval^ A & , had to be decreased as R 
increased to keep the process convergent. For R^ very small, the 
iteration interval had to be made very small again to get meaningful 
result (proper divergence from the zero deflection line). 
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To eliminate much of this difficulty, equations 39 and 40 were 
developed. A similar iteration procedure to that of equations 47 was 
initially used to investigate the new equations, and it v/as found that 
tm iteration interval of 0. 01 radians gave satisfactory results for 
between 0. 01 and 100. 

Eqtaations 39 and 40 appeared to be the more satisfactory of 
the two sets of equations for programming on a digital computer, since 
the starting solutions of numerical methods are closely related to the 
simple iteration procedure. 

A Fortran statement for the solution of equations 39 and 40 
was v/ritten utilizing a subroutine prepared by Mr. K. H. Redner of 
Computer Sciences Corporation. The subroutine applies the method 
of Runge-Kutta- Gills as a starting solution, and the Adams-Moulton 
predictor - corrector formulas as a continuing solution. The facilities 
of the Computer Center at the California Institute of Teclinology were 
used in writing the program, and the IBM 7090 computer at the Jet 
Propulsion Liaboratory v/as used for the acUial solution. A Fortran 
statement is included as Appendix A. 

To investigate the general effect of changes in "S' and n , 
Stokes drag law v/as assumed to apply, and solutions were obtained 
for six values of between 0. 01 and 50 for three sets of assumed 

constants. 

1) Mo- J~Z , Y=h4 ^ n - 2 .o , 

2) Mo ='/F' , y = /.4 , n= 1.0 , 

3 ) Mo-yr , r= F2 , n^ 2 .o ■ 
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The e>daaust conditions of a typical large rocket were assumed 
as initial conditions for the investigation of the effect of allowing 
to vary according to equation 23. The initial conditions chosen were; 



jX 0 



Po 



Go 

V 

n 



3, O , 
9.5 5 X 10 



-7 I bf - Sec 






O. 0)38 



Ibnn 



3,100 ^ps (Tq~ 2820“R , 
\.Z5 , 

/. O . 



Solutions were obtained for of 0. 1 , 1.0, and 50 for two par- 

ticle sizes, 0“ =1 micron and cj" = 2. 0 microns. 
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V. RESULTS AND DISCUSSION 

The resxilts of this investigation are plotted on Figs. 3 through 
10 for 0° to 70° of flow deflection angle. Flow deflection angle, S , 
was chosen as a more meaningful parameter than -O' for plotting and 
analyzing the results since it is determined solely by the physical 
boundaries of the flow. The fan angle and deflection angle can be re- 
lated by the usual Prandtl-Meyer functions. 

Figure 2 is a graphical representation of equation 6 giving Ro 
as a function of and 0" for particles having a density of 

210 Ib/ft (as for ). Initial gas characteristics plotted are for 

air at 100°F and 1000°F and for the assumed rocket exliaust conditions 
given in the previous section. is not greatly ciffected by a change 

in initial gas conditions, but is primarily a fxmction of and g" 

For particles between 1. 0 (a and 4. 0 u in diameter and initially dis- 
placed by several cm. from the flow boundary (as in the case of a 
rocket exhaust), Ro will normally be greater than 10. 

The data plotted in Figs. 3, 4, and 5 is based on the assump- 
tions that Mo ~ v/^ , y =• /. 4- , n ^ 2 . O , and Stokes drag law ap- 
plies with tlie factor -f - O . por these conditions, particles with 
Ro 5 O will have little deviation from the streamline, and slip 
velocities will be low. For Ro ^ 0.0 1 there will be large tan- 
gential slip and very little deviation from its original flow direction 
through the first 25° of gas deflection. Radial slip is small compared 
to tangential slip, being a maximum of . 164 for — O- Ol at about 

25° flow deflection. Near this point the maximum relative Mach 
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number is . 993 and masiimum Reynolds number for a 4. 0 4 particle 
is 8 . 38. B'or - j o, tbe majiimum radial slip is less tlian one per 
cent while the maximum Mach number is only 0. 151 and maximum 
Reynolds number for 4. 0 u particle is 2. 02. The Mach numbers and 
Reynolds number botii drop sharply on either side of the maximum 
values. 



The assumed viscosity - temperature function, 

n 



t. 



il) 






equation 19 

(19) 



is a well-established approximation that applies over a limited tem- 
perature range when n is assigned the constauit value appropriate to 
this range. By allowing n to take its expected limiting Vcilues, its 
effect on the solution can be observed. FigTirss 6 and 7 show curves 
for which Stokes drag, I. 4- , and ■= \f~^ v/here assumed for 

n - I. o and n = 2. o . A small cheinge in trajectory path and tan- 
gential slip velocity and a negligible cheinge in radial slip velocity 
were observed. Thus, little error is introduced by the selection of 
an average value of n over an extended temperature range. 

A change in "ZT from 1. 4 to 1.2 does have a fairly large ef- 
fect on the trajectory path v/hich is due principally to the change in 
the gas streamline. The tangential slip velocity is Increased also, 
which tends to increase the path deviation from the streamline, but 
this is a small effect. Radial slip velocity is not appreciably affected. 

The drag correction factor, -f , v/hich allows a departure from 
Stokes drag lav/ to compensate for slip flow or free molsciilar flow 
effects, was proposed by Rannie (ref. 3 ) as an analogue to a similar 
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development for heat transfer coefficients. Although it has not been 

verified experimentally, it does give proper results for the extreme 

values of i. o and ^ y. o » shovJd give a good 

Re Re 



approximation for intermediate values. The factor 



Ms 



can be 



shown to be (ref. 7 ) 



=: o. 4 9 9, 



S 



Tt y 2<r 



where is the mean free rRolec^^lar path of the gas. Thus, if 

^ >> | . O , free molecular flow e:dsts, and if << /. o 

f?e ^ ^ 

continuum flow exists. The a.rea bebvveen these limits is the slip 

flow regime. 

From equation 27 we have, for \- Z and n - o , 



^ ^ f COS c< \IQ 

\,Cos(X9 + °^)/ 



(48) 



For air or rocket exliaust gases at moderate temperatures, and for 
particles the order of a micron, k will lie between 0. 03 and 0. 3 . 
Thus, a particle of this general size, passing through a Prandtl- 
Meyer e>q>ansion, will start in the slip flow regime and pass into free 
molecvilar flow at large erqjansion angles as the enclosed cosine fionc- 
tion in equation 48 becomes large compared to unity. 

Figures 8, 9, and 10 show the effect of ^ on tlie trajectory 
patlis and slip velocities for 1. 0 micron and 4. 0 micron diameter 
j>articles. The deviation from the uncorrected Stokes drag law tra- 
jectories is large for particles as small as 4. 0 micron diameter. 
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amd increases markedly as the particle size is reduced to 1. 0 micron. 
In the case of a 1. 0 micron diameter pax-ticle, an = 50 does not 
closely follow the streamline as it did for Stokes drag, and Ro - o. \ 
now shows little deflection from its original flow direction. For 
Rq = O. 1 a maximum relative Mach number of 0. 924 and ma>dmum 
Reynolds number of 6. 73 is reached for a 4. 0 micron particle. This 
compares to maxlmxims of 0. 853 and 6. 5i respectively for Stokes 
drag. 

The maximum values of Reynolds nuirbers obtained are greater 
than the limits set on equation 23 for the determination of -f . A 
Reyirolds number of 8. 0 would change tlie equation to 

^ = O. 65 + 4.0 5 vT — ’ 

T Re 

(ref. 3 ), which would make a 23 per cent error in -f for a Mach num- 
ber of 0. 95 . While this is an appreciable error, it is effective over 
a small portion of the total expansion and diminishes rapidly on both 
sides of maidmum , so it sho^^id not greatly affect the overall 



results. 



• CONCLUDING REr./IARKS 



vl 



In the flow regime where Stokes drag law applies, the tra- 
jectories of spherical particles passing through a Prandtl-Meyer ex- 
pansion fan can be specified by a single dimensionless parameter, , 
wliich is a function of initial gas conditions, initial radius, and particle 
diameter. If the Pveynolds number is greater than about 0. 5, tiae tra- 
jectories are also a function of Re . V.'hen the mean free patli of the 
gas molecules is of the same order or large compared to the particle 
diameter, Stokes drag law must be modified to correct for tlie discon- 
tinuous nature of the flow, and the trajectories then become functions 
of three parameters, , Re , and 2vi^/Re . 

It v/as found that the temperature dependence of viscosity need 
not be precisely known to ensure satisfactory results, but changes in 
the specific heat ratio, y , do appreciably :-ffect the trajectory paths. 

The radicd velocity lags are in general very small compared 
witli the tangential velocity lags, and can, in many cases, be neglected 
completely. 
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APPENDIX A 

FORTRAN STATEMENT 



CI^ENSICK Y(IC)» YCCTO), C(10), C(IO) 

CIKENSICN Y(10), YCCTO), C(IC), C(IC) 

5C0 FCP^AT (7F10.A) 

51C FCR^'AT (ICHl GANKA = E13.6, 13H MZERC = E0.6, 

19H N - F5.2, 13H RZERC = E13.6, 13H DELTA = Eli. 6//) 

IC READ INPUT TAPE 5,5CO, CGAPPA, CMZERO, CN, RZEKD,DELT, C^'RE, CO 
KRITE CUTPUT TAPE 6, 510, CGAPPA, CHZERO, CN, RZERO, DELT 
Y(l) = l.C 
Y(2) = 1.0 
Y(3) = O.C 
T = 0.0 
KAY = I 

CLAMCA = SCRTFC (CGAPMA-1.0 )/(CCAPPA+1.0) ) 
e = n .0/ (CLAKCA#*2 ) )-CN+l .0 

ARC = SORTF (1.0♦CLA^^DA*CLAPDA*(CMZERG•CMZERC-1.0) ) 

ALPHA = ACCSF( l.C/ARG) 

C(l) = ( l.C/CLAPCA)-CLAPDA 

CC2) = C( I )/CLAPCA 
C(3) = CCSF (ALPHA) 

C(A) = CLAKCA*(C(3)**(B-1.0) )*RZER0 
C(5) = 2.U/CLAPDA 
C(6) = CMRE / CC 

tiRITE CUTPUT TAPE 6 , 5C0 ,CL AKCA , B , ARG , ALPHA 

VsRITE CUTPUT TAPE 6 , 5 CC ,C ( 1 ) , C ( 2 ) , C ( 3 ) , C ( A ) , C ( 5 ) , C(6) 

CALL CEC(K,3,T,Y,YDCT,.Cl,5.CE-7) 

GC TC (lC0,20C,2CO,30C),K 
ICC C(l) = CLAHCA*T+ALPHA 
C(2) TANF(C(D) 

CO) = CCSF(C( 1) ) 

C(A) = ( l.C/ (C( 3) **B) ) 

C(5) = l.C/ ( l.0+C(6) • ((C(3)/C(3)) **(B-l.O))) 

YCCT(l) = Y(3) 

VCCT(3) = (C(5)) »Y(3) «D(2) - C ( 2 ) * Y ( 1 ) *D ( 2 ) *D ( 2 ) 
l-(C(A)»Y(3)/( (Y(l)#*3)*Y(2)*Y(2)) )*D(A) * D(5) 

YDCT(2) = “C ( I ) ^YC 2) *C (2 )+C( A ) #D( A ) *( I .0/ ( ( Y( 1) »*3) »Y(2 ) ) 

1-1 .0/ ( Y ( 1 ) *Y( I ) ) ) *C(5 ) 

IF (KAY-2) 150, 160, 16C 
160 call CEQ2 

15C V.RITE CUTPUT TAPE 6, 500, C(l)i D(2), C(3), D(A),D(5) 

WRITE CUTPUT TAPE 6, 60C 

6CC FCRMAT (//6H THETA, 9H X, lAH Y, 

llAH Z, 16H V/VS,IAH U/US, 

21AH R/RC, I3h MS, 13H RE//) 

KAY = 2 
GC TC 160 

2CC Y(A ) = Y( 1 )*Y (2) 

Y(5) = Y(A)*(1.0-CLAMDA*YDCT(1)/(Y(1)*D(2))) 

Y(6) = ( (C(3)/C(3) )**(C(2)+1.0) )/Y( 1) 

Y(9) ( ( l.C-Y( A ) ) »*2 )+ ( ( ( l.C-Y( 5 ) )*D (2 )/CLAMDA )**2 ) 

Y(7) = SCRTF(Y(9) ) 

Y(8) = Y( 7) *( (C( 3) /CO) ) ** (B-l.O) ) 

WRITE CUTPUT TAPE 6, 61C, T, ( Y ( JACK ), JACK = 1,B) 

610 format (/F7.3,8EIA.6) 

2C5 IFU-3. 0)210, 210,220 
21C CALL CEQl 
GC TC 205 
220 GC TC 10 

CIMENSICN Y(IC), YDCT(3), C(10), D(10) 

3CC WRITE CUTPUT TAPE 6, 62C 

620 FCRMAI (22H ERROR RETURN FRCM CEO) 

CC TC 10 

ENC( 1,0,0,0,0,0,C,0,0,1,0,C,0,C,0) 
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